﻿using System.Diagnostics;
using System;
using System.Windows;
using System.Windows.Controls;
namespace IMap
{
    public static class Core
    {
        #region Beijing_1954 CS
        //Geographic Coordinate System:

        //Name: GCS_Beijing_1954

        //Alias:

        //Abbreviation:

        //Remarks:

        //Angular Unit: Degree (0.017453292519943299)

        //Prime Meridian: Greenwich (0.000000000000000000)

        //Datum: D_Beijing_1954

        //Spheroid: Krasovsky_1940

        //Semimajor Axis: 6378245.000000000000000000

        //Semiminor Axis: 6356863.018773047300000000

        //Inverse Flattening: 298.300000000000010000
        #endregion

        //Clarke 1866 Spheroid
        private const double SemimajorAxis = 6378206.4;
        private const double SemiminorAxis = 6356583.8;
        private const double InverseFlattening = 294.9786982;
        private const double MinLatitude = -85.05112878;
        private const double MaxLatitude = 85.05112878;
        private const double MinLongitude = -180;
        private const double MaxLongitude = 180;
        private const double PI = Math.PI;
        private const int TILE_SIZE = 256;
        private const double AVERAGE_R = 6370996.81;//地球平均半径 

        public static double NumClip(double n, double minValue, double maxValue)
        {
            return Math.Min(Math.Max(n, minValue), maxValue);
        }
        public static int NumClip(int n, int minValue, int maxValue)
        {
            return Math.Min(Math.Max(n, minValue), maxValue);
        }

        public static PixelPoint LatLong2Pixel(LatLong latlong, int level)
        {
            return Mercator2Pixel(LatLong2Mercator(latlong), level);
        }

        public static LatLong Pixel2LatLong(PixelPoint pixelPoint, int level)
        {
            return Mercator2LatLong(Pixel2Mercator(pixelPoint, level));
        }

        public static MercatorPoint LatLong2Mercator(LatLong latlong)
        {
            double latitude, longitude;
            latitude = NumClip(latlong.Latitude, MinLatitude, MaxLatitude);
            longitude = NumClip(latlong.Longitude, MinLongitude, MaxLongitude);
            //http://blog.csdn.net/xcl119xcl/article/details/5791821

            double K = SemimajorAxis;
            double f = 1 / InverseFlattening;
            double e = Math.Pow((2 - f) * f, 0.5);

            double B = latitude * PI / 180;
            double L = longitude * PI / 180;

            double S = (1 - e * Math.Sin(B)) / (1 + e * Math.Sin(B));
            double mercatorX = K * L;
            double mercatorY = K * Math.Log((Math.Tan(PI / 4 + B / 2) * (Math.Pow(S, e / 2))));

            Debug.Print(@"墨卡托投影坐标(X=" + mercatorX.ToString() + @" Y=" + mercatorY.ToString() + @")");

            return new MercatorPoint(mercatorX, mercatorY);
        }

        public static LatLong Mercator2LatLong(MercatorPoint mercatorXY)
        {
            double mercatorX, mercatorY;
            mercatorX = mercatorXY.X;
            mercatorY = mercatorXY.Y;

            double K = SemimajorAxis;
            double f = 1 / InverseFlattening;
            double e = Math.Pow((2 - f) * f, 0.5);
            int iterativeValue = 10;
            double B = 0;
            for (int i = 0; i < iterativeValue; i++)
            {
                B = 2 * Math.Atan((Math.Exp(mercatorY / K)) * Math.Pow((1 + e * Math.Sin(B)) / (1 - e * Math.Sin(B)), e / 2)) - PI / 2;
            }

            double L = mercatorX / K;
            double latitude = B * 180 / PI;
            double longitude = L * 180 / PI;
            return new LatLong(latitude, longitude);
        }

        public static PixelPoint Mercator2Pixel(MercatorPoint mercatorXY, int level)
        {
            double mercatorX = mercatorXY.X;
            double mercatorY = mercatorXY.Y;
            double pixelX = (mercatorX / GroundResolution(level));
            double pixelY = (mercatorY / GroundResolution(level));

            Debug.Print(@"像素坐标（X=" + pixelX.ToString() + @" ,Y=" + pixelY.ToString() + @")");

            return new PixelPoint(pixelX, pixelY);
        }

        public static MercatorPoint Pixel2Mercator(PixelPoint pixelPoint, int level)
        {
            double pixelX, pixelY;
            pixelX = pixelPoint.X;
            pixelY = pixelPoint.Y;
            double mercatorX = pixelX * GroundResolution(level);
            double mercatorY = pixelY * GroundResolution(level);
            Debug.Print(@"墨卡托投影坐标(X=" + mercatorX.ToString() + @" Y=" + mercatorY.ToString() + @")");
            return new MercatorPoint(mercatorX, mercatorY);
        }

        public static Tile Pixel2Matrix(PixelPoint pixelPoint)
        {
            int matrixX, matrixY;
            matrixX = (int)(pixelPoint.X / TILE_SIZE);
            matrixY = (int)(pixelPoint.Y / TILE_SIZE);
            Tile matrixXY = new Tile(matrixX, matrixY);
            Debug.Print(@"瓦片编号[" + matrixX.ToString() + @"," + matrixY.ToString() + @"]");
            return matrixXY;
        }

        // Determines the ground resolution (in meters per pixel) at a specified latitude and level of detail.
        public static double GroundResolution(int level)
        {
            return Math.Pow(2, 18 - level);
        }

        public static double Distance(LatLong latlon1, LatLong latlon2)
        {

            //R*arccos[cosβ1cosβ2cos（α1-α2）+sinβ1sinβ2];                          
            double beta1 = (latlon1.Latitude) * PI / 180;
            double beta2 = (latlon2.Latitude) * PI / 180;
            double alaph1 = (latlon1.Longitude) * PI / 180;
            double alaph2 = (latlon2.Longitude) * PI / 180;
            double dist = AVERAGE_R * Math.Acos(Math.Cos(beta1) * Math.Cos(beta2) * Math.Cos(alaph1 - alaph2) + Math.Sin(beta1) * Math.Sin(beta2));
            return dist;
            //地图上距离，而非球面距离
            //double dist = Math.Sqrt(Math.Pow(myline.X1-myline.X2,2)+Math.Pow(myline.Y1-myline.Y2,2))*Core.GroundResolution(this.Level)/1000;            
        }

    }

    public class IImage : Image
    {
        const int TILE_SIZE = 256;
        public IImage()
        {
            this.Height = TILE_SIZE;
            this.Width = TILE_SIZE;
            this.Stretch = System.Windows.Media.Stretch.UniformToFill;
            // this.Stretch = System.Windows.Media.Stretch.None;这样设置有问题
            //由于部分图片为72dpi，例如u=x=1650;y=444;z=13;v=014;type=web&fm=44.png
        }
    }

    public struct MercatorPoint
    {
        public double X, Y;

        public MercatorPoint(double x, double y)
        {
            this.X = x;
            this.Y = y;
        }

        public static MercatorPoint operator +(MercatorPoint mpt1, MercatorPoint mpt2)
        {
            return new MercatorPoint(mpt1.X + mpt2.X, mpt1.Y + mpt2.Y);

        }

        public static MercatorPoint operator -(MercatorPoint mpt1, MercatorPoint mpt2)
        {
            return new MercatorPoint(mpt1.X - mpt2.X, mpt1.Y - mpt2.Y);

        }

        public static MercatorPoint operator *(MercatorPoint mpt, double scale)
        {
            return new MercatorPoint(mpt.X * scale, mpt.Y * scale);

        }

        public override string ToString()
        {
            return "MercatorPoint (X=" + X.ToString() + ", Y=" + Y.ToString() + ")";
        }

    }

    public struct PixelPoint
    {
        public double X, Y;

        public PixelPoint(Point pt)
        {
            this.X = pt.X;
            this.Y = pt.Y;
        }
        public PixelPoint(double x, double y)
        {
            this.X = x;
            this.Y = y;
        }

        public static PixelPoint operator +(PixelPoint ppt1, PixelPoint ppt2)
        {
            return new PixelPoint(ppt1.X + ppt2.X, ppt1.Y + ppt2.Y);

        }

        public static PixelPoint operator -(PixelPoint ppt1, PixelPoint ppt2)
        {
            return new PixelPoint(ppt1.X - ppt2.X, ppt1.Y - ppt2.Y);

        }

        public static PixelPoint operator *(PixelPoint ppt, double scale)
        {
            return new PixelPoint(ppt.X * scale, ppt.Y * scale);

        }

        public override string ToString()
        {
            return "PixelPoint (X=" + X.ToString() + ", Y=" + Y.ToString() + ")";
        }


    }

    public struct LatLong
    {
        private double latitude;
        private double longitude;

        public double Longitude
        {
            get { return longitude; }
            set { longitude = value; }
        }

        public double Latitude
        {
            get { return latitude; }
            set { latitude = value; }
        }

        public LatLong(double latitude, double longitude)
        {
            this.latitude = latitude;
            this.longitude = longitude;
        }

        public static LatLong operator +(LatLong ll1, LatLong ll2)
        {
            return new LatLong(ll1.Latitude + ll2.Latitude, ll1.Longitude + ll2.Longitude);

        }

        public static LatLong operator -(LatLong ll1, LatLong ll2)
        {
            return new LatLong(ll1.Latitude - ll2.Latitude, ll1.Longitude - ll2.Longitude);

        }

        public static LatLong operator *(LatLong ppt, double scale)
        {
            return new LatLong(ppt.Latitude * scale, ppt.Longitude * scale);

        }

        public override string ToString()
        {
            return "LatLong (Latitude=" + latitude.ToString() + ", Longitude=" + longitude.ToString() + ")";
        }

    }

    public class Tile
    {
        public double PixelX, PixelY;
        public int Column, Row;

        public Tile(int column, int row)
        {
            this.Column = column;
            this.Row = row;
        }

        public Tile(double pixelX, double pixelY)
        {
            this.PixelX = pixelX;
            this.PixelY = pixelY;
        }

    }

    public struct TileMatrix
    {
        internal int MinX, MinY, MaxX, MaxY;
        internal int Width, Height;

        public TileMatrix(double width, double height, PixelPoint pixelPoint)
        {
            BBox bbox = new BBox(width, height, pixelPoint);
            Tile lefttopMatrix = Core.Pixel2Matrix(bbox.lefttopPoint);
            Tile rightbottomMatrix = Core.Pixel2Matrix(bbox.rightbottomPoint);
            this.MinX = lefttopMatrix.Column;
            this.MinY = rightbottomMatrix.Row;
            this.MaxX = rightbottomMatrix.Column;
            this.MaxY = lefttopMatrix.Row;
            this.Width = MaxX - MinX + 1;
            this.Height = MaxY - MinY + 1;
        }
    }

    public struct BBox
    {
        internal double width, height, minX, minY, maxX, maxY;
        internal PixelPoint lefttopPoint, leftbottomPoint, righttopPoint, rightbottomPoint;

        public BBox(double width, double height, PixelPoint centerPixelPoint)
        {
            this.width = width;
            this.height = height;
            minX = centerPixelPoint.X - width / 2;
            minY = centerPixelPoint.Y - height / 2;
            maxX = centerPixelPoint.X + width / 2;
            maxY = centerPixelPoint.Y + height / 2;
            lefttopPoint = new PixelPoint(minX, maxY);
            leftbottomPoint = new PixelPoint(minX, minY);
            righttopPoint = new PixelPoint(maxX, maxY);
            rightbottomPoint = new PixelPoint(maxX, minY);
        }
    }

    //Baidu的算法墨卡托投影转换为经纬度
    public class BaiduMC2LL
    {
        //http://hi.baidu.com/feiyyu/item/448a59485809e2a960d7b936
        double[] MCBAND = new double[] { 12890594.86, 8362377.87, 5591021, 3481989.83, 1678043.12, 0 };

        double[,] MC2LL = new double[,]{ 
                                    {0.000000014105261721162549, 0.00000898305509648872, -1.9939833816331,200.98243831067961, -187.2403703815547, 91.6087516669843,-23.38765649603339, 2.57121317296198, -0.03801003308653, 17337981.2},
                                    {-0.0000000074358563895655369, 0.0000089830550977262389, -0.78625201886289,96.326875997598464, -1.85204757529826, -59.369359054858769,47.400335492967372, -16.50741931063887, 2.28786674699375, 10260144.86},
                                    {-0.000000030308834608988258, 0.00000898305509983578, 0.30071316287616,59.742936184422767, 7.357984074871, -25.383710026647449,13.45380521110908, -3.29883767235584, 0.32710905363475, 6856817.37},
                                    {-0.000000019819813049305521, 0.0000089830550997795349, 0.03278182852591,40.316785277057441, 0.65659298677277, -4.44255534477492,0.85341911805263, 0.12923347998204, -0.04625736007561, 4482777.06},
                                    {0.00000000309191371068437, 0.0000089830550968121549, 0.00006995724062,23.109343041449009, -0.00023663490511, -0.6321817810242,-0.00663494467273, 0.03430082397953, -0.00466043876332, 2555164.4},
                                    {0.000000002890871144776878, 0.0000089830550958054071, -0.00000003068298,7.47137025468032, -0.00000353937994, -0.02145144861037,-0.00001234426596, 0.00010322952773, -0.00000323890364, 826088.5}
        };

        public Point ConvertBaiduMC2LL(Point cE)
        {
            double[] cH = new double[MC2LL.GetLength(1)];
            Point cF = new Point(Math.Abs(cE.X), Math.Abs(cE.Y));
            for (int i = 0; i < MCBAND.Length; i++)
            {
                if (cF.Y >= MCBAND[i])
                {
                    for (int k = 0; k < MC2LL.GetLength(1); k++)
                    {
                        cH[k] = MC2LL[i, k];
                    }
                }
            }
            Point t = Convertor(cE, cH);
            return t;
        }

        private Point Convertor(Point cF, double[] cG)
        {
            double T = cG[0] + cG[1] * Math.Abs(cF.X);
            double cE = Math.Abs(cF.Y) / cG[9];
            Double cH = cG[2] + cG[3] * cE + cG[4] * cE * cE + cG[5] * cE * cE * cE + cG[6] * cE * cE * cE * cE + cG[7] * cE * cE * cE * cE * cE + cG[8] * cE * cE * cE * cE * cE * cE;
            if (cF.X < 0)
            {
                T *= -1;
            }
            else
            {
                T *= 1;
            }
            if (cF.Y < 0)
            {
                cH *= -1;
            }
            else
            {
                cH *= 1;
            }
            return new Point(T, cH);
        }
    }

}

